/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of LINAL (LINear ALgebra)
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/
#include "gpu_test.hxx"

int main(int argc, char **argv) {

  linal::Flat_    A,  C, D, norm;
  linal::Hier_   hA, hC;

  FLA_Obj *alpha, *beta;

  if (argc != 10) {
    printf("Try :: %s ", argv[0]);
    printf("[n_thread] [n_thread_gpu] [n_gpu] ");
    printf("[n_bin] [n_work] ");
    printf("[uplo] [trans] [alpha] [beta]\n");
    return -1;
  }

  int verbose = 0;

  int
    n_thread     = atoi( (argv[1]) ),
    n_thread_gpu = atoi( (argv[2]) ),
    n_gpu        = atoi( (argv[3]) ),
    n_item_bin   = atoi( (argv[4]) ),
    n_item_work  = atoi( (argv[5]) );

  int uplo, trans;
  int datatype = TEST_DATATYPE, n_device;

  cudaGetDeviceCount( &n_device );
  // ---------------------------------------
  // ** Initialization   
  FLA_Init();

  if ( atoi( (argv[2]) ) ) uplo  = FLA_UPPER_TRIANGULAR;
  else                     uplo  = FLA_LOWER_TRIANGULAR;
  if ( atoi( (argv[3]) ) ) trans = FLA_TRANSPOSE;
  else                     trans = FLA_NO_TRANSPOSE;
  if ( atoi( (argv[4]) ) ) alpha = &FLA_ONE;
  else                     alpha = &FLA_MINUS_ONE;
  if ( atoi( (argv[5]) ) ) beta  = &FLA_ONE;
  else                     beta  = &FLA_MINUS_ONE;

  printf("** TEST ENVIRONMENT **\n");
  printf("N         = %d\n", N);
  printf("GPU_MN    = %d\n", GPU_BMN);
  printf("CPU_MN    = %d\n", CPU_BMN);
  printf("N thread  = %d\n", n_thread);
  printf("N device  = %d\n", n_device);
  printf("N gpu     = %d\n", n_gpu);

  // ---------------------------------------
  // ** Matrices
  A.create   (datatype, N, N);
  C.create   (datatype, N, N);
  D.create   (datatype, N, N);
  norm.create(datatype, 1,1);

  FLA_Random_matrix(~A);
  FLA_Random_matrix(~D);

  hA.create(A, GPU_BMN, GPU_BMN);
  hC.create(D, GPU_BMN, GPU_BMN);

  FLA_Copy( ~D, ~C);

  if (verbose) {
    FLA_Obj_show("=A=", ~A, "%4.3f", "===");
    FLA_Obj_show("=C=", ~C, "%4.3f", "===");
    FLA_Obj_show("=D=", ~D, "%4.3f", "===");
  }

  // ---------------------------------------
  // ** FLAME
  FLA_Syrk( uplo, trans, *alpha,
            ~A, *beta, ~C);

  // ---------------------------------------
  // ** GPU
  omp_set_num_threads( n_thread );
  linal::init_gpu(n_thread_gpu, n_gpu, n_item_bin, n_item_work );
  linal::set_cpu_blocksize(CPU_BMN);
  linal::set_work_balance(1);
  linal::set_computing_model( TEST_COMPUTING );
  linal::create_flat_gpu( datatype, GPU_BMN, GPU_BMN );

#pragma omp parallel
  {
#pragma omp single nowait
    {
#pragma omp task 
      linal::dense::syrk( uplo, trans, *alpha,
                          hA, *beta, hC);
    }
  }

  linal::disp_gpu_data_access(stdout);
  linal::disp_operation_stat(stdout);
  
  linal::set_computing_model( LINAL_CPU );
  linal::free_flat_gpu();
  linal::finalize_gpu();

  if (verbose) {
    FLA_Obj_show("=~C=", ~C, "%4.3f", "===");
    FLA_Obj_show("=~D=", ~D, "%4.3f", "===");
  }

  // ---------------------------------------
  // ** Check
  FLA_Triangularize(uplo, FLA_NONUNIT_DIAG, ~C);
  FLA_Triangularize(uplo, FLA_NONUNIT_DIAG, ~D);

  FLA_Axpy(FLA_MINUS_ONE, ~C, ~D);
  FLA_Norm1( ~D, ~norm);

  if (verbose)
    FLA_Obj_show("=Diff D=", ~D, "%4.3f", "===");


  float  norm_f = *((float *)FLA_Obj_buffer_at_view(~norm));
  double norm_d = *((double*)FLA_Obj_buffer_at_view(~norm));

  int rval;

  printf("- TEST::");
  for (int i=0;i<argc;++i)
    printf(" %s ", argv[i] );
  printf("\n");

  switch (datatype) {
  case LINAL_SINGLE_REAL:
    if (norm_f < LINAL_ERROR_TOL) {
      printf("PASS::Norm :: %E \n", norm_f);   rval = 0;
    } else {
      printf("FAIL::Norm :: %E \n", norm_f);   rval = -1;
    }
    break;
  case LINAL_DOUBLE_REAL:
    if (norm_d < LINAL_ERROR_TOL) {
      printf("PASS::Norm :: %E \n", norm_d);   rval = 0;
    } else {
      printf("FAIL::Norm :: %E \n", norm_d);   rval = -1;
    }
    break;
  }

  // ---------------------------------------
  // ** Matrices
  A.free(); hA.free();
  C.free(); hC.free();
  D.free();

  norm.free();

  // ---------------------------------------
  // ** Finalization
  FLA_Finalize();
  return rval;
}



